Nagaoka ferromagnetism in the two-dimensional infinite-^7 Hubbard model 
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We present different numerical calculations based on variational quantum Monte Carlo simulations 
supporting a ferromagnetic ground-state for finite and small hole densities in the two-dimensional 
infinite-?/ Hubbard model. Moreover, by studying the energies of different total spin sectors, these 
calculations strongly suggest that the paramagnetic phase is unstable against a phase with a partial 
polarization for large hole densities 5 ~ 0.40 with evidence for a second-order transition to the 
paramagnetic large doping phase. 
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The Hubbard model was originally introduced to ex- 
plain ferromagnetism in strongly correlated transition 
metals 

H = -t cl^Cj^a + UY^ rii^tn^^i, (1) 

where ( ) stands for nearest neighbors, ci^a (cjg.) de- 
stroys (creates) an electron with spin a at site i, and 
n-i^a = cJ^Ci^cr- In the following all energies are measured 
in units of t. In this simple and appealing Hamiltonian 
the nontrivial ingredient is that pairs of electrons with 
opposite spins that occupy the same site pay an on-site 
energy U. For a large Coulomb interaction U, electrons 
can prefer to occupy only a single spin sector, up or down, 
thereby minimizing the strong Coulomb repulsion. In- 
deed, as Nagaoka showed in a milestone paper j2|, a sin- 
gle hole in the infinite-?/ Hubbard model on any finite 
bipartite cluster with periodic boundary conditions, for 
any dimension d > 2 has a fully polarized ferromagnetic 
ground-state, i.e., the ground-state has maximum spin 
S. Unfortunately, the Nagaoka theorem refers to a sin- 
gle point in the phase diagram, which is not thermody- 
namically significant. Therefore, a key issue in strongly 
correlated electron systems is to understand if a ferro- 
magnetic ground-state can survive at finite hole densi- 
ties, finite Coulomb repulsions, or can be stabilized by 
slight modifications to the pure Hubbard model. 

Recently there is a renewed interest in this field H] 
and many authors have focused on the two-dimensional 
infinite-C/ Hubbard model but the situation is still 

controversial. In particular, Putikka et al. using the 
high temperature expansion, found that the ground-state 
does not have a saturated magnetization for all hole den- 
sities. Due to the uncertainties of their extrapolation 
at zero temperature, their calculation cannot be consid- 
ered a definite proof of the complete absence of a Na- 
gaoka ground-state. On the other hand, on any finite 
lattice with L sites, the total spin S of the ground-state 
is strongly dependent upon the number of holes and the 
boundary conditions chosen. For instance for two holes 
on any finite cluster with periodic boundary conditions 
it is possible to show 0] that the ground-state spin is 



not maximum, whereas for certain finite number of holes 
such that the ferromagnetic ground-state is nondegener- 
ate with periodic boundary conditions, the ground-state 
has maximum spin [p|jic|]. In general, if the boundary 
conditions are such that the ferromagnetic (paramag- 
netic) state is a nondegenerate closed shell, ferromag- 
netism (paramagnetism) is energetically preferred. Given 
the above difficulties, we believe that a reasonable ap- 
proach to study the possible ferromagnetic instability in 
the infinite-t/ Hubbard model is to study a finite density 
of holes and periodic boundary conditions such that the 
singlet state is nondegenerate, thus favoring the singlet 
state. Whenever, for large size systems, we find that a 
state with a finite spin polarization has a lower energy 
than the singlet one, despite the fact that in the for- 
mer case the closed shell condition may be not satisfied, 
this clearly implies a ferromagnetic instability weakly de- 
pendent on boundary conditions. In order to proceed 
with this scheme, large cluster sizes are necessary since in 
small systems there are few and irrelevant singlet closed 
shell dopings. In the present work we are able to con- 
sider fairly large lattice sizes so that our findings rely 
mainly on closed shell dopings, even though some of the 
calculations were performed with open shells, in order to 
interpolate between two closed shell densities. 

The Quantum Monte Carlo method (QMC) has re- 
cently been a subject of intense activity |ll|-[l^ due 
to its ability to handle the exponentially large Hilbert 
space dimension of strongly correlated systems and it is 
one of the few reasonable approaches for large size sys- 
tems. In this letter we have made extensive use of a 
recent QMC technique [Q, which is particularly suit- 
able to estimating the ground-state energy of a Hamil- 
tonian H by applying a few Lanczos steps to a given 
starting variational wave function Iv^g)- By using the 
Green-function Monte Carlo with stochastic reconfigu- 
ration it is possible to compute the best energy state 
I = Sfc=o (^kH^ I , which defines the wave function 
with p Lanczos steps applied to \^g)- This can be done, 
even for large size, at least up to p = 2. Expectation val- 
ues of the Hamiltonian {'^p\H\'^p) and corresponding en- 
ergy variance cr^ = [{^p\H'^\^p) - (*p|i7|^'p)2] /L'^, can 
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then be computed by QMC, with high statistical accu- 
racy. As well known, an accurate estimate of the energy 
per site can be obtained by computing the energy as a 
function of the variance for different wave functions [ p^ . 
The variance is zero for an exact eigenstate and therefore, 
if different variational wave functions are systematically 
approaching the ground-state, the exact value of the en- 
ergy can be estimated by extrapolating the different cal- 
culations up to the zero variance limit. The systematic 
approach of the Lanczos technique to the exact solution 
is remarkably accurate even for large system sizes. Thus 
the variance extrapolation technique combined with the 
Lanczos algorithm represents, at the moment, one of the 
most simple and efficient schemes to estimate the ex- 
act energy or, at least, the error of our best variational 
energy. As shown in Ref. iQ, with a reasonably good 
variational state, it is possible to estimate exact energies 
within the statistical error, even for L ^ 100. 
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FIG. 1. (a): variational energies for p = 0,1,2 in the 
singlet subspace for 28 electrons on a one-dimensional chain 
of 50 sites. The continuous line is a quadratic fit. The ex- 
trapolated point with zero variance is also shown, (b): vari- 
ational energies for p = 0, 1,2,3 in the singlet subspace for 
10 electrons on a two-dimensional cluster of 18 sites. The 
continuous (long-dashed) line is a linear (quadratic) fit of the 
p = 1,2 (p = 0,1,2) results. The extrapolated points with 
zero variance and the FNLS results are also shown. 

For the present study, the Lanczos method also has the 
remarkable advantage that, if \^g) has definite quantum 
numbers, the Lanczos method will conserve them exactly, 
and therefore it is possible to compute the lowest energy 
in the different spin sectors. In order to show the quality 
of our best variational wave function on large size sys- 
tems, we will also use the standard fixed-node approxi- 
mation ]T6[ |, by using the "p — \ variational state as the 
guiding wave function, which will be denoted by FNLS. 
This approximation typically leads to slightly lower ener- 
gies than the best p = 2 Lanczos ones. However, within 
the fixed-node approach, the total spin is not conserved 



and therefore this approximation is not useful for esti- 
mating the spin dependence of the energy, which is our 
main task. Moreover, within the fixed-node approach, 
it is not possible to calculate the energy variance and 
therefore to assess the accuracy of the calculation. 
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FIG. 2. Differences between the singlet energy Es=o/L 
and the fully polarized one Efm/L for the 98-site lattice as 
a function of the hole density for different approximations: 
p = 1 (open triangles) and p = 2 (full squares). The re- 
sults for the FNLS approximation in the Sz = subspace 
(open circles) and the extrapolated values (full circles) are 
also reported. The lines are guides to the eye. In the inset: 
Es=o I L for the 50-site lattice (full circles) and for the 98-site 
lattice (empty squares) are reported. The ferromagnetic en- 
ergy for the 50-site lattice (dashed line) and for the 98-site 
lattice (continuous line) are also shown. 

The p = variational wave function is 

|*g) = VnVgJU (l + /fe4,T^-fca) (2) 

where Vn is the projector onto the subspace of N parti- 
cles, Vg is the Gutzwiller projector, which forbids dou- 
bly occupied sites, J = expij j ''^ij^i^j) i^ ^ J^^" 
trow factor, defined in term of the hole density at site i 
hi = (1 — ?T.i|)(l — rii^i), and 7 is a variational parameter. 
For the potential v, we take the exact analytic form deter- 
mined by considering the holes as hard-core bosons at the 
same density JTrf . The variational parameters fk of the 
pairing wave function may also represent the Gutzwiller 
wave function in the particular case of /j, = 0(eF — ea:), 
where e^^ is the free Fermi energy, ek is the correspond- 
ing dispersion ek — ^2t{cos + cos ky) and Q is the step 
function. For closed shell fillings, this wave function is a 
singlet. For open shell, we take f^. with a small d-wave or 
s-wave symmetry in order to split the bare degeneracy, 
the resulting wave function being a singlet. Analogously, 
by using the particle-hole transformation on down spins, 
we are able to consider a number of up electrons which is 
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different from the number of down electrons, correspond- 
ing to a finite total spin S > 0. We use the 45° tilted 
squares with L = l\f2 x sites with periodic boundary 
conditions, in order to have the full spatial symmetries 
of the infinite lattice and remove the huge degeneracy of 
the conventional I x I square lattices near half-filling. 

For the infinite-C/ Hubbard model, we have found that 
the best variational wave function of the class given in 
Eq. (||) has a free electron determinant, /fe = 9(6^ — e/c), 
even for nonzero total spin. The Jastrow term J provides 
a sizable improvement in the variational energy (of order 
0.01 per site). Definitely, this class of wave functions is 
particularly accurate in the large doping limit where the 
paramagnetic state is expected to be stable apart from 
Kohn-Luttinger superconducting instabilities, which are 
however expected to lead to negligible energy gain and 
are irrelevant for our purposes. In the more delicate 
small-doping region, based on variational QMC calcula- 
tions fisll , it appears that if nontrivial order parameters, 
other than the ferromagnetic one, show up in the ground- 
state, they should lead to a marginal energy gain. Within 
this assumption, even in this case, the Gutzwiller wave 
function can be considered a good starting point for our 
projection technique. 

In order to test the Lanczos step technique and apply 
it to the infinite-C/ Hubbard model, we show in Fig. ^the 
energy as a function of the variance for two closed shell 
cases where the exact values are known: a small two- 
dimensional 18-site lattice, and a larger one-dimensional 
chain. The approach to the exact ground-state energy is 
remarkably good, even when, as shown in Fig. |l|(a), the 
p = variational wave function is about 20% off in energy 
in the one-dimensional case, where the exact solution 
is known by Bethe ansatz. For small two-dimensional 
systems, the p = variational wave function is zero 
for several electronic configurations of the Hilbert space, 
whereas, as soon as we apply the first Lanczos step, this 
pathology is substantially removed. Thus, in this case, 
it is natural to fit linearly only the p = 1, 2 results, the 
extraplated value being not affected by the inclusion of 
the p = 3 point (which is possible for this small cluster), 
see Fig. |l|(b). Nevertheless, even by fitting quadratically 
the p = 0, 1, 2 points, we obtain a rather accurate result. 
On larger size systems, L ^ 50, we have verified that 
the number of zeros in \^g) is negligible and that all the 
three energies can be safely included in the same fit. 

We perform a systematic study of the spin instability 
and we consider the 50-site and the 98-site lattices. The 
results for these two systems are both in qualitative and 
quantitative agreement, suggesting that, for L ~ 100, the 
finite size corrections are small compared to the energy 
difference between the singlet and the saturated ferro- 
magnet. In Fig. |2| we report the difference in energy be- 
tween the singlet and the fully polarized state for different 
hole densities and QMC techniques for the 98 sites. For 
large doping the homogeneous singlet state is energeti- 



cally well below the ferromagnet but, by decreasing the 
number of holes, there is a clear evidence of a fully polar- 
ized ground-state. Indeed all the Monte Carlo methods 
(variational with p = 0, 1, 2 and FNLS) predict that the 
singlet energy is below the fully polarized one only for 
(5 > ^SF, where the critical doping 8sf ~ 0.28-;- 0.30. Of 
course, 8sf gives an upper bound for the critical density 
(5ci , below which the saturated ferromagnetism is stable. 
Indeed, unsaturated ferromagnetism can appear below 

^SF- 
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FIG. 3. Variational energies with p = 0, 1,2 in different 
total spin sectors for 42 electrons on 50 sites: 5* = (trian- 
gles, continuous line), S = % (squares, long-dashed line), and 
S = 16 (circles, short-dashed line) . The fully polarized energy 
for the same lattice is also reported. 

If we consider the extrapolated value of the energy, 
6sp is pushed at much lower doping ~ 0.20. This 
transition point appears to be in very good agreement 
with the best variational estimate for the instability of 
the fully polarized state |^,^. In our numerical calcula- 
tion we are able to study directly the instability of the 
singlet-paramagnetic state, starting from the "safe" large 
doping limit. The fact that we have obtained a value 
for the transition point very similar to the one found in 
Refs. |6[0], supports the validity of both approaches in 
predicting a stable ferromagnetic region. 

In order to show that the fully polarized state has the 
lowest energy in the low-doping regime, in Fig. ^ we show 
the variational energies for 42 electrons on 50 sites, cor- 
responding to (5 = 0.16 < 8sF- We consider all the pos- 
sible values of the spin which fulfill the closed shell con- 
dition, whereas the fully polarized state is not a closed 
shell. This condition, as discussed in the introduction, 
certainly frustrates the fully polarized state and favors 
the singlet. Moreover, we have considered a linear fit of 
the p = 0,1,2 energy results, yielding a slightly lower 
extrapolated energy compared to the one obtained with 
the more accurate quadratic fit. Even with this choice, 
the ferromagnetic state remains clearly stable. Moreover, 
the energy decreases monotonically with increasing total 
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spin. The outcome of this calculation is that, as we in- 
crease the system size, the ferromagnetic state becomes 
less and less frustrated. For the singlet state instead, the 
boundary conditions do not play a crucial role and the 
corresponding energy per site is very weakly size depen- 
dent (see inset of Fig. ^). 




FIG. 4. (a): variational results with p = (full trian- 
gles), p = 1 (full squares) and p = 2 (full circles) for the 
ground-state energy as a function of the total spin S for 58 
electrons on 98 sites. The extrapolated energies (empty cir- 
cles) are also shown. The ferromagnetic energy is marked by 
the arrow and the lines are guides to the eye. (b): the same 
for 66 electrons on 98 sites. 

Since in the low-doping region the fully polarized state 
is found to be stable, it is reasonable to expect a partially 
polarized ground-state for intermediate doping. For very 
large doping, 6 > 0.50, the ground-state is found to be 
paramagnetic and the energy is monotonically increasing 
with increasing total spin, implying a finite spin suscep- 
tibility. For every variational Monte Carlo calculation 
the results are consistent with a paramagnetic ground- 
state and a finite spin susceptibility. By decreasing the 
hole density, the low-spin states collapse to the singlet 
ground-state and eventually become degenerate with it. 
For 40 holes on 98 sites, we found that, for all the Monte 
Carlo techniques used, the singlet is degenerate, within 
the statistical error, with the 5 = 8 state, both the states 
being closed shell, see Fig. ^(a). Indeed, although in the 
variational calculation there is a small but significant dif- 
ference between the energy of the S — and the 5 = 8 
state, by increasing the accuracy, with p Lanczos steps, 
the two energies tend to become closer and closer and 
even the extrapolated values give a degeneracy of these 
two spin sectors. This indicates that at doping ~ 0.40 
the uniform spin susceptibility diverges and a finite mag- 
netic moment sets in for S < 6c2- These results strongly 
support a second-order transition to a partially polarized 
ferromagnet, and a nonzero magnetic moment is in fact 
seen at doping S ~ 0.33, where the energy is minimum 



for 5 - 10, see Fig. |(b). 

To conclude, unsaturated ferromagnetism begins al- 
ready at ^ ~ 0.40 and rather reliable numerical evidence 
is given that Nagaoka ferromagnetism is not simply an 
irrelevant property of a single hole, but is realistic phys- 
ical behavior of strongly correlated electrons. Since at 
small but finite doping ferromagnetism is stabilized by a 
macroscopic energy gain, this property is robust against 
small perturbation of the model, such as a large repulsion 
U OT a small positive coupling J in the t — J model . 
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Preliminary results confirm that for 5 ~ 0.1 and J/t < 
0.05, the ground-state is ferromagnetic, whereas, by de- 
creasing 5, first ferromagnetism melts, and then, at very 
small doping, antiferromagnetism is finally stabilized. 
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